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ABSTRACT 


Master's Thesis which discusses nature of allocation 
problems. Danskin Algorithm for solution of a convex func- 
tion to be minimized sewer a closed conwex set is developed. 
An example of application involving solution of a 3600 
variable allocation problem using a computer is provided. 
Paper includes analysis of solution and discussion of pro- 


blems encountered in Computer application... 
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I. INTRODUCTION 


A large selection of problems faced within the mili- 
tary and commercial environments consist of finding the 
"best" way of using available resources. In mast cases 
the "best" method of using these resources is that which 
mEeovides the most desirous return to the using agency or 
meoanazation. Typically, this return cam be represented 
mera payotf associated with the use of resources, or in 
mathematical terms it is a function - the value of which 
is dependent on and determined by the amount or value of 
resources expended. Naturally, the availability and use 
of these resources are bounded and, therefore, act as a 
Gonstraint on the using organization. This constraint 
Gauild be that the organization is required to use some set 
faninmum amount of a particular resource, or in the more 
usual sense only a given maximum anount of the resource 
is available for use. In the former case the resource is 
termed lower bounded, and in the latter case the resource 
is considered upper bounded. In many such problems the 
resource is both upper and lower bounded and usually some 
allocation of the resource between the bounds is acceptable. 

If the situation exists where the bounds on a resource 
are known but where varying amounts of the resource can be 
expended in a number of different ways, then the manner in 
which the resource is allocated becomes of critical importance 


womenic ali@catine organization im terms of the payort se 
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receives from choosing a particular allocation scheme. 

By simple use of a Crystal ball. zandom selection device. 
or some other means; and by exercising care not to violate 
the bounds, the organization can adequately allocate. How- 
ever, it cannot be sure it has found the best allocation - 
Seneme . 

If the payoff realized can be expressed ag a mathe- 
fetical function of the weseurces allocated, then the pro- 
blem becomes susceptible to solution by mathematical 
programming methods. In many cases then the "best". or 
optimal allocation of resources can be found and proven to 
be optimal. These type problems can be described as "al- 
iecation type problems. Ttypmedlly, the valte representa- 
tion of a resource variable is a non-negative number less 
@ian one and indicating the percent of the total resource 
available which is to be expended in the particular manner 
@r activity associated with the given variable. 

ime formulation of Whe payoff functions usmally the 
meecidine factor in determining what type of mathematical 
programming method could~be used to find the optimal al- 
mocatiren or optimal solution to the problem. For instance, 
Mmamtie payott function is@¥inear in nature, then 1t may be 
possible tomsolve the probilenmin USiIn?g exrstine Minear pre. 
gramming methods. However, many payoff functions cannot be 
expressed as linear relationships even when attempts are 
made at reasonable approximations. In this case one must 


turn to other methods. In some cases it may be true that 








no known mathematical method for optimizing CxS tS.) -uneetnac 
case it may be proper to reformulate the problem or possibly 
Cecurna tO theseryseolebaae 

This paper will address a particular category of al- 
location problems, develop a method for solving problems 
which fit into this category, and provide a detailed exam- 
mle OL a practical problem solved by the me tnodedce veloped. 
eredit for development of this method belongs te Dr. John 
M. Danskin, Professor - U. S. Naval Postgraduate Scheol. 
This method, employing an aleori chin termed the Direction 
Finding Algorithm, was presented by Dr. Danskin in various 
milassroom Veetures presented during his Gerure at that 


eenool . 


A. PRESEN wRnON Ob SHE SRO EM 

Consider an allocation type preblem in Which the pay- 
Ser function is convex for concave) and it isedecided that 
the optimal solution is that solution which minimizes (or 
Maximizes) the payoff function. Now, consider possible al- 
location of resources to be represented as X = (X1,X2,-.-X)) 
mere X; = percent of resources expended in activity #1, 
m= percent of resources expended in activity #2, ete.; for 
me- |... ,N. 

As indicated earlier, certain constraints must be heeded 
in finding the optimal solution. These constraints can be 
formulated as follows: 

(1) It is necessary and desirable to use up all of the 


Besoupees available, 








«a 


(2) It may or may mot be mecessary tomisetsene ner 
cent of the resource in any given activity, 

(3) The percent of resource expended in activity i 
cannot exceed its upper or lower bounds, and 

(4) Negative use of the resource is meaningless and 
not possible. 

This problem as described can then be mathematically 


expressed as: 


maximize f(x) where f(x) is concave 
Or 


minimize f(x) where f(x) is convex 


subject to constraints: 


where a; and b. BeMreSent wie Upper and lowerepeundse Ge - 
eeectively of the allocation to the jth eG te ee teny oe reemmte tel 


follows that the following relationships are true: 


Essential to the development and use of the Direction 
Finding Algorithm is the fact that the payoff or objective 
function must be differentiable. Therefore, this require- 
Ment that the gradient €x1st6 and can be calculated Ws eon. 


cIdened eas an added Constraint. 








This Category of sallocation problems is = sees came dase 
the category for whieh the aboue- to pe develencdalarocte non 


Finding Algorithm can be used to find the optimal solution. 





TI. DERIVATION AND DEVELOPMENT OF THE ALGORITHM 


In this section the Direction Finding Algorithm will 
be developed and a step-by-step procedure Heer aned fe Onte 
mes application, 

Recall the original problem is to maximize (minimize) 


a concave (convex) differentiable function of the vector 


X = (X1,X25-++5X)) with constraints as follows: 
ie Xeo= 1 0 <a.<b. 
: il - ji 1 
a 
de Xe eepe eeecee <i < Obie 
1 ee ‘i 1 , i 


The basis of the algorithm is to find the direction 
iemtastest increase. This méans it is mecessary to maxi- 
mize the Directional Derivative which is defined in Ref. [1] 


as 


MCS) 22 se 8) (Ge 
u i i i 


where 2 represents the gradient and y represents the direc- 
tion of the gradient. 
To develop the algorithm, the following constraints are 


miposed on y, the direction vector: 


eae” 
1 
yY¥% = 1 
een] 
1 

Y; 2 Paes xX, = a; 


y. < O 1f x. = b. 





The entire problem then becomes 


maximize ZR*¥y 


subject to I Y; = ] 


2 Nee = 0 
Yo 0 if X; = a, 
ves O if xX, = bs 


moth the assumption that the maximum 1s@men—-negacive. 

Mme first step to solwame this preblem will be to show 
Mecessary conditions for solution by applicationget the 
well-known Kuhn-Tucker (K-T) Conditions which can be found 
imienet, [2] and many other places in recent literature on 
non-linear programming. The constraints can be rewritten 


and LaGrange Multipliers assigned as follows: 


CONSTRAINT MUI PLIER 
7 1 
2 Ys td A 
= 2 11 
ii yy; < 0 d 
? 
i : Pi 
~ ' ? 
eM Ss 8 ae 
SR (0) Scr eee Tor 
SL 1 1 i 
ve < 0 tt Xe = De 1 
a a a al 


Application of the Kuhn-Tucker Theorem then results in 


V (Q*°y) = A'V(Z ve =) + A''V(1-Zy=) 
SUNG y J (ey) 
ae) er ica Gian) 


10 





where 


V aS usual representc te jar acsretine 


LE’ amplies Summation Over tiesewe selon wien 
ee = a 
1 1 
X'' implies summation over those i's for which 
x. =' Dee 
1 1 


If at this point only the ith component of the gradient 


is considered, the above equation leads to 


. = ¥i(2y! + oi! Sy! =a ay 
oA reer Bt ) ial Pi wa F 
where again Vv, enters consideration only if xX, = a. and 7; 
is considered only if x, = b,. The Ys represents the ith 


Gomponent of the y vector in the optimal solution as pro- 
Vided for in the Kuhn-Tucker Theorem. Letting A = 2\A' - 2A"! 
and yp = u' - u'', the following results can readily be de- 


auced from the previous equation and the original con- 


Seraints: 
0 e- 
r ee if a, Sass < b. 
Or Itaee= 4-0 and youu 
1 il il 
Q. =< ee a 
‘i or if Xs b. and ee 0 
0 . = On 
r Y; +4 V5 if Xs as and 4 0 
0 : _ 0. 
r Y; +H + 7, if x, b. and Ys 0 


SInNGe ve 2 0 and T, 2 0 from the stim —-luiewer lheoremr 


mes cqQuivalentily truce that 


11 
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0 7 = 0 _ 
Qs < A Y; +4 if xX. b. and Ve 0 (1) 
0 . = 0 ol 
Q. cae Veil ia: Xx. a. and Ya 0 (29 
= ene 1) otherwise conditions satis- 


Dene the original con- (3) 
SEraAalmes. 

Therefore, if the optimal solution exists, then A and 
uw eXist and are non-negative and satisfy the above condi- 
enon Ss. 

Me mext step 1s £0 €onsider the suffirereneaict sence 
conditions. It is asserted that if X} and uw exist and i is 
mem-necsative and if y° Satisties results in (1), G2). and 
(3) above, then y’ maximizes IQ+y. 

To prove this let y be any arbitrary direction vector 
Satisfying the conditions derived via the Kuhn-Tucker 
Theorem. Observe that Q-y = IM.y;, = (2. - u)y, since 
LY = = 0 from the original problem constraints. This is true 
for all i representing elements of the 2 and y vectors. 
Therefore, the summation can be broken into three groups 


fae cerms aS indicated in (1), (2). and (3) above; i.e., 


wey = Et (Q,-u)y, : Ett (Q.-w)y, + phtT(Q.-w)y; 
where ea and ptt see SUItat el ONSmeMe: Ceiiis. Tie) (co) 
and (3) ie ope Celina a 0 eo bnee (2. - UL) < 0 from the K-T 
mecessary condition (1) and Yee 0 from the original con- 


meraint. By similan seasoning It Can be sccm that eo Se 


WMiemeioGe. 1 telS enue ethat wey en 


vill 0 


However, rT (Q.-u)y, = } Yaa from K-t cond erem 


1 
eneae ‘ a ae 
fo) and AL Yi, =A all Y;¥; since y; = Gor all aor 


= 12 








ligigi » 0 yl o2 eee 
contained in 2% », and A 44 Ys; § Mail " rat a eatin x 


1 1 ii 
from the Schwartz inequality and from the original con- 


straints. Therefore, Q*y < A, but Qey° = E(Qe-n) ye = ALys* 
=. And since X > 0, it is true that Q-y < Q+y° with the 
resultant fact that y’ maximizes 2*y as was to be shown. 
Furthermore, the maximum value is A or in terms af the suf- 
ibeLlency Conditions; if X can be found then the cptimal 
solutit1on 15 @btained. 

Ae thisepoimt Items comvemient to rewrite eamaeron 


(3) and the conditions necessary for the equality: | 


= 0 1 
2. AY = + U ie a. < Xs < bs 
Omeeet Sond Crandon 
i qi at 
Om Xie are ee Oe 
i i ie 


Letting 2' indicate a summation over the terms ful- 


filling these conditions and summing, results in 


! = 2,0 t a ! 
> 2. AX Y; toe ul ae 
oT (4) 
a ee 
y= N! z 2 


mmere N' is the total number of terms fulfilling the condi- 
tions. If the equation is squared first and then summed, 


me Expression 
hos 16s sy)? (5) 


mesmits. [t also follows from substitution that 
0 _ 7 
V5 (25 ey (6) 
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With a method of calculating values for nu, A, and Y; 
now in hand, development of the algorithm can proceed. 

For notational ease, let DS = the "Distinguished Set" 
and be defined as that set of indices i such that a.<x,<b., 


or X. = a. and Ys > DOO X 2 


= b. and y? < 0. In other 
al i i i i 


words, the distinguished set represents those elements of 
the allocation vector, X, where the element does not lie on 
the boundary; or if the element is on the boundary, then 
the corresponding direction vector:element implies a move 
away from the boundary. With this notation equation (4) 


@on be rewritten as 


Ws a ieps % (7) 

where Nyc is the total number of elements in DS... 

the elements Contained in DS can be separated into 
mearee Categories as follows: 

Category #1. The X element is not on its boundary and 
tmere is no restriction on the y element. 

‘Category #2. The X element is on its lower boundary 
and the y element is positive. From equation (6) then 
a. oe Vie 

Category #3. The X element is on its upper boundary 
and the y element is negative. From equation (6) then 
5 eT < 

It is evident now that the algorithm for finding y°, 


or the direction of fastest increase, must search out and 


find each of the elements included in the three categories 
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above. Recall that uw 1s actually seme wo ei ace Or eal 25 One 
which Y; # 0. This means it is necessary to determine u 
concurrently with determining the ellieinei Gs aimee) >. neete meno tee 


owing procedure prevades the method) tor dome just nase 


A. THE DIRECTION FINDING ALGORITHM 


ceed z Z 
Step #1. Let uw = y— Gino %, and DS = cons tart. 
DS 
Steep 72, Place into 8S each index deeer ewe a,<xi<b,. 
Nete that this implies no restriction on Y;- Note also that 
throughout the algorithm, once i has been entered into DS, 
imcmne longer considered aSeamcandidate Lommenutim.: 


Step #5. af DS = 1¢}, then temporarily let wy= mint, }. 
| i 


Do not place this © into DS (the minimum gradient eFement 
serves simply as a starting point). If DS # {o}, then cal- 
mikate uw from equation (7). 

DkLep #4. Scamen indices fon cases where os bs. Find 
the smallest 2. associated with this search. If this 25 <wiie, 
@iemeenter this i imto DS and recalculate w. If this 2 2 4, 
Maem proceed to Step #6. 

Sieepet 5. Repeatmatep #4 until no waddatienal 1°s can 
me entered into DS. 

Seer O. ocarch Iitivees monmertses where xX; = a. Find 
the largest 2 ASSOCIated withmunts scaren.s If this ae ea 
then enter this i into DS and Toca] cuneiie i cee ele stele 2. lees 
then proceed to Step #8. 

Scope 7-7 @kepeadl stem FOrunti | Normadds tromeima, sameciml 


be entered into DS. 
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Step #8. Repeateste pt 4) Sei io eons ee eee eee mecicle 
ditional i's can be entered into DS. 

Step #9. Calculate A = [ 2% (2. - nu) 27%. If >} = 0 

1e Ds 

aon LOD 

Step #10. Calculate y° from equation (6). 

This concludes the algorithm. At this point two 
Special case considerations should be discussed: 

CASE I. The possibility exists that only one element 
will be entered into the DS. In this event = 0 and the 


optimal solution lies on the boundary. y’ 


at. tis pone 
indicates that any better solution would violate the orig- 
inal constraints. 

CASE™11. “lieimere than oneeienene has been entered 
into DS and A = 0, then 2. = uw for all i. In ids eVeme 
all 2. are obviously equal and the optimal solution has 
been reached. In fact, the value of A is a measure of how 
elose the procedure has come to the optimal solution. It 
fo be true in a practical application of the algorithm that 
X will diminish more slowly with each successive computation 
as the optimal solution in neared. It may be necessary in 
this case for the user to make a determination of the pre- 
emslon desired and to terminate calculations when a satis- 
factory degree of precision has been reached. 

In the next section an example application of the algo- 


rithm in a practically oriented problem is provided. 
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Til. EXAMPLE APPLICATION OF DIRECTION FINDENG ALGOR SIEM 


Consider a military oriented problem in which an ord- 
nance delivering organization has responsibility for and 
the capability of firing its ordnance into a particular 
sector of territory. Assume that this organzzation has 
some means of assigning probabilities to the potential 
existence of a target at any of a number of specific loca- 
tions within its sector of responsibility. Assume also 
that the organization knows what damage it can expect from 
a specific type round delivered at a specific point. Know- 
m2 this information, theserganization would then like to 
Few how to optimally distribute the inventory it intends 
to deliver so that the probability of target survival is 
minimized.! 

WSe a £IySt Step togs@lution of this premiiem, the  sceucn 
of responsibility is assumed to be square in geometric shape 
and an imaginary grid is imposed upon the sector. The in- 
dividual squares of the grid are numbered sequentially be- 
ginning in the top left corner and working first across and 
then down. Location of a target (or the possible location 
@f a target) can then be referenced by a specific numbered 
square of thewerid- Expectucd —daRace Inman), soda SCUa re seam 


be calculated from the number of rounds which impact in that 


This problem is formulated and presented for solution 
Mimatnests Paper prepared by Captain William Ay Hessen, 
UsiiGew sateen. S. Naval Postgraduate School in Pebruary, 
March 1971. 


Ne 





square and/or in nearby squares. ~Allocaticnuso. soOndnanec 
is represented by the number of rounds delivered to a 
Srceri1¢e rrdesquanre. 
whe sumvival pmobabagiey, of gagitanme tm. ar the objective 
function which is to be minimized, can now be expressed as 
a mathematical equation 
EGY) : — : Ee oe 


where the variables have interpretations as follows: 


Pp, = probability that target is located in ith square, 
Bs = probability Of destruction of a tamget appeanane 
in the jth Square by sa round détenmating ian the 


jth SqGtlare (this as reterred to as the damage 
ine ea OTe. 

Ver = Percent Or tOtLal inventory tombemcemenacd alien 
cated to the jth square, 

Fay) 


survival probability of the target and is a 
LUncCtOnmon themalilocatiron Of it. 

Note Jue tise pOlmltnthat the Survival = Une «momen (1 — 
i a convex function and the variable Y is actually a 
closed convex set (there are only a finite number of rounds 
which can be fired). Also note that in the trivial case 
momen Y = (00. ...,0), t-e., no firing takes place, the 
Survive time ClOmencaices. to F (i)ie= Fé De = I which 1s ap- 
propriate . 


The problem now is to determine what allocation of Y 


will minimize F(Y). 
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It can be seen at this point “that amm@order to obtarn 
realistic and usable results , WG iacmiie cess a Ope eco. ec! 
with a large number of squares. Correspondingly, the allo- 
eatton vectoreris composedsot a large mimmbeieon elements. 

The problem can now be represented as follows: 


minimize F(Y) = ; me 


SUbJeCt CO VeonSeuaunts 


2 = 1.0 
1 


and 


From previous discussion and assuming the Ps and oe 
aoe Known beforehand, this problem can be solved using the 
ierecti0on Finding Algorithm. Observet@Wat the S@rst partial 
@erivative required in the use of the algorithm is 

FY) = ; P; Ba a See ee 

SOluULHenM Of eha sep vob lem Vida themDirec hronmsic ice Looe 
rithm was accomplished on the IBM 360 computer system at the 
U. S. Naval Postgraduate School. Programming was done in 
FORTRAN IV Language using the G-Level Compiler. A copy 


of the program listing is included after the appendices. 


A. EXPLANATION OF EXAMPLE 

SInGcesoIespUrpese Of thisesexample ms) to demons ueacemaa. 
plication of the algorithm and not to solve a given existing 
PEOD Lem ww asstgenment of values to parameters or chne propmcem 
Was based on academic interest and convenience rather Chan 


Peteetediily “orlented . 
19 





A grid size of 3600 squares, or equivalently a 60X60 
matrix was used. It is felt this grid size selection is 
consistent with what might be used in a truly practical 
application and at the same time it enlarges the problem 
sufficiently so that information about machine time for 
calculation would be meaningful. 

The distribution of probabilities (p; "s) used in this 
example was uniform in nature. A p = 1/3600 was assigned 
to each of the 3600 grid squares. ‘In a normal application 
it is expected that most grid squares would have ap = 0.0 
assigned to them and only a relatively small number would 
be assigned a positive probability. Those positive assign- 
ments would be most likely based on observer reports, known 
routes of movements, specific areas offering good cover and 
concealment, and other considerations such as these. How- 
6ver, in this example a uniform distribution was selected 
wor two reasons. Firstly, making all p's positive 1s”"con- 
Sidgered a worst case condition with respect to the number 
of individual calculations required and the resultant ma- 
chine time necessary for execution. Since solution by the 
Direction Finding Algorithm is an iterative process, machine 
time is an important consideration and a known upper limit 
on ae factor is extremely meaningful. Secondly, a uniform 
distribution was used so that effects on the boundaries of 
the grid could be observed. Since the damage function (P55 3) 
implies interaction between grid squares, it is not immedi- 
ately obvious what the optimal allocation along the edges of 


the grid should be. 
20 
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Selection of an appropriate damage function would be 
based on the type of weapon being employed and the result 
of statistical data obtained for that weapon. Inherently, 
the damage function selected directly affects total machine 
mame necessary ta solve the problem. “Howevetmeenis is a 
variable which is dependent on a specific situation and 
mee OL prime consideration in this paper. A simple expon- 
eieial relationship was selected for the examplesbecause 
ae wasc tele tiat an Exponential function would be = aepre-— 
meiitaeive Of most damage functions in terms of accuracy 
and machine time necessary for calculation. 

Given the above explanation of Drobleneeamameters Le 
now becomes necessary to develop a computer oriented logi- 
cal approach to solution of the problem. 

ive Voctealecavence=tolfowed an) finde ene @diree tron 
of fastest increase is explained earlier in this paper. At 
this point it should be noted that the program actually em- 
ploys two Direction Finding Algorithms. [In one case 
(SUBROUTINE DF$ALL in the program listing) the y or direc- 
@mome element associated with each and éveryselement of the 
allocation vector is calculated. This 15 necessary = petone 
a final optimal solution can be obtained. However, while 
Pork ine toward the optimal solution, a significant reduction 
mm Calculation time can begrealized by finding the direction 
of fastest increase only on the variables which have a posi- 
tive allocation currently associated with them or a positive 


Vet case the aiocation 1s currently zeroveewor thismreason 


Zi 








the program maintains a list (SUBROUTINE LIST) of those ele- 
ment indices which meet this criteria. With each tteration 
then, SUBROUTINE DF$LST is employed to find the direction 
of fastest increase only for those elements appearing on 

the "list."' When the optimal solution is found for the 
euerent “list, execution of the programeshifewmeeescUb— 
ROUTINE DF$ALL. After execution of DF$ALL, a new list is 
formed and the above outlined process is repeated. Only 
when optimization is indicated through the use of DFSALL 
iseeXeCCUbLIon terminated. 

In order to employ the Direction Finding Algorithm as 
indicated, various other calculations obviously must be 
Made. Other routines in the program are written to accom- 
plish tasks such as evaluating the gradient, evaluating 
mac ODJECctive tumetion, and revising the current aivocaeren 
to one which is closer to the optimal solution. These tasks 
must necessarily be carried out during each iteration of 
the solution process. A general flow chart outlining the 
logical procedure followed is depicted in Appendix A. A 
complete listing of subroutines employed along with a des- 
Sription of the purpose of each is provided in Appenddeesar 

Resultsmebtained inethe solution to this cape prob- 
lem ane included. It 1S interesting to note in these results 
that the optimal solution has no allocation of resources in 
the rows or columns immediately adjacent to the edges of the 
grid. Obviously the damage resulting in these grid squares 


from rounds allocated to and impacting in nearby interior 


Le 
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grid squares is considered suffticiemeusremeemp noate for loscmo. 
damage effect that would occur if rounds were allocated to 
the outermost rows and/or columns. The higher allocation in 
those nearby interior grid squares suggests this conclusion. 

It is also interesting to note that the optimal solu- 
tion is not completely symmetric as might be expected from 
using a uniform type distribution of the p;'S.. Inspection 
Of the final allocation shows the upper left corner of the 
array somewhat different from allocations assigned in the 
meoer threevcorners: tihalseis the result of ‘themstartsng )al- 
location used in the solution process and the precision con- 
sidered acceptable within the program. The initial starting 
feasible solution was for the entire inventory to be assigned 
to the extreme upper left grid square. Th®s heavy allocation 
i just one location prevented a lesser posative ailocation 
from appearing in nearby grid squares as the program stepped 
along towards optimality. The program was actually moving 
toward a totally symmetric final allocation, but settled for 
something less because of the precision cutoff programmed 
into the problem. If a higher degree of precision were de- 
manded and if the program were rewritten to provide greater 
machine accuracy, the program would tend to a completely 
ametric OVATE Hues so) oGe alone 

Certain difficulties were encountered in programming 
this problem for machine solution. A discussion of some of 
these difficulties and comments on observations is considered 


Seeacademme Interest to the reader at ethis point. 
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Sinee use of thewmirection Finging Galeer tenes an 
iterative process which moves from any feasible starting 
solution toa better ene and then? wheriaee tor an approxi 
mate optimal solution, selection of an initial starting 
wont Can be a Critical matter. “initralwatcempes at solving 
mrs problem were dome by assigning the ~mtrre inventory 
bor tne first grid Square and proceeding from there. Anal- 
ysis of results obtained from this "start from scratch" ap- 
proach indicated that machine time*in exeess of twelve hours 
would be required to completely solve the problem. . This re- 
quired time was considered unacceptable and, therefore, an 
alternate approach was devised. The sector covered with a 
3600 grid was initially considered to be covered with a 400 
element grid. An optimal solution was fotmd for the 400 grid 
Case and this solution was then used as a starting solution 
for the 3600 grid case. Results in terms of machine time 
were remarkable. The optimal solution to the 400 case was 
mieeained in 17.5 seconds of machine execution time. Over- 
all machine time to the final optimal solution was six 
minutes nine seconds. 

fepolme Of Meerest=imethe Program aiseuhne value of "'D" 
eetne distance which the allocatiom is®maved itm each itera- 
eon . ‘Recall the lower and upper boundaries of the alloca- 
tion elements are 0.0 and 1.0 respectively. The lower bound 
on D is 0.0 which is synonomous with no move at all. The 
upper bound on D is the maximum distance across the Y-space 


and is ¥ 2. In some calculations within the program D is 
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controlled by a variable which is near its boundary and 
trying to Move in he eda cet wom so fame eon date elu 

result of this situation in large problems such as this 
example is that D becomes very small and this leads to 

short moves on each iteration. When many variables are 
ellose to their botndaries and trying to move sin trac dinec- 
pron. the overall ampact is a signiticamt mMmig@rease in 
Moehine execution time. However, there are various places 
in the solution process where D can be successfully reset 
Eoua predetermined value so that time 19 mOt wasted on 
relatively small moves when larger moves are permitted. 
There is also the problem of what to do with D when move- 
Ment 1S made in a direction of increase but it has not re- 
sulted 1n am wmproved value of thie objeetmwe function. To 
mesolve these difficulties, an initial D value of 0.5 was 
selected. This had the effect of making the first step a 
Benge one. A reset value of 0.1 was decided upon and each 
mime anew list was compaled, D was reset to thas value. 

A modification procedure was used to adjust D when a move 
was made but an improvement in the objective function was 
met realized (rightly enough referred to as a "failure" ): 

In this case D was eae cut in half and thesmove tried 
again. in the case of Success in improvang the value of the 
objective function, D was left alone and the same value used 
himEhe Next iteration Ehewonlysexcepetomemomunis procedure 
was in the case of a "success" which had been immediately 


Preceded Dyas taditure. in this Circumstance I was halved 
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for the next iteration. To understand the rationale be- 
hind this consider the following simplified two-dimensional 
representation of a concave function which is to be maxi- 


mized: 


we tepoilnt A represent the current allocation and the wdiree— 
tion of fastest increase is obviously to the right. Sup- 
pose the present D is equal in value to (B-A). When the 
move to B is attempted, a “failure” o¢€eurs since the value 
of tite function is less than it was at p@int A. The pro- 
eeaure now calilis for halving D and try to move™again. This 
mates the move is to point C and a “success” is realized. 

mem 1S Not adjusted at this point, the mext iteration 
would call for an attempted move back to point A which was 
Ehe starting point for the previous iteration. To preciude 
@ais return to a previous allocation, D is cut in half 
after the occurrence of a success which had been immedi- 
ately preceded by a failure. It should be noted that this 
procedure becomes particularly useful when the optimal 
Solution is being neared and oscillations around the optimal 


SOlUuemon Start OCCuring 
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Another consideration worth sotneGgem1s themint) uence 
of machine accuracy on@the results of eeieeproblem. Sina 
problem of this sone with potentially geo evar iables pes - 
tive simultaneously, the effects of calculating, rounding, 
and cumulatively adding these variables must be taken into 
account. Early attempts at solution were made using single 
precision under FORTRAN IV. This approach was quickly 
enanged to double precision for all real vamrables used in 
the program. Even so, results obtained were considered 
suspect in accuracy beyond iia T@ tellustrate “thempmob- 
Mem, in this example the inpatral P; value assignment to 
each grid square was made by assigning each grid square a 
value of 1/400. This value is supposed to be accurate in 
FORTRAN IV to ie As a check device the program then 
sums all p's assigned and prints the results. The re- 
sultant sum of 400 separate additions differs from 1.0 by 
2.24x10 °. Whether or not accuraCyeat thaselevelis mod 
erable is up to the user, but most certainly is worthy of 
some consideration. 

Machine storage space and execution time requirements 
in this problem were sufficiently significant to be worthy 
of comment. Because of the requirement to maintain a large 
number of large arrays of stored information (in the common 
meorage region alone thewerare sixX @mmays each of which con- 
Tinos OO0Ndolib le precistenyreal numbers) thissexampllienpno- 
gram run on the IBM 360 System required 352,000 bytes of 


Storage. Some space saving techniques were employed and, 
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undoubtedly additionalimenes could @be tolmd. Gilewever, it 
should be recognized that the storage requirements are not 
a direct result of the methodology used in the Direction 
Finding Algorithm, but rather the result of the nature of 
mae problem being solved. Asementionedsprevrousis, total 
execution time necessary was Six minutes nine seconds. The 
approach of solving the problem on a small scale first and 
then using these results as a starting solution for the 
mare er Case was an instrumental factor in keeping execution 
time within reason. Other techniques involving the applica- 
tion of "clean coding" contributed toward reducing time 
requirements. One technique developed is considered to be 
of academic interest and involved the calculation of geo- 
metrical distance between squares. This calculation is 
necessary in evaluation of the gradient and in evaluating 
the function value. The distance calculation involves sim- 
my measuring the herizontal and vertical distance between 
the two grid squares in question and then evaluating the 
‘alanre root Of the sum of squares. To illustrate the Signa- 
ficance of this simple routine, it is done in excess of one 
foeellion times im each iteration of the program. The execu- 
muon time Gequired to aati the built-in machine square root 
function and evaluate the argument is known to be in excess 
of 200 microseconds. Since this obviously contributes to 
exorbitant machine execution time, it was decided to make 
the distance computation between any two grid squares once 
and for all at the beginning of the program and then store 


this information for later reference as required. This was 
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possible since the distance evaluations were somewhat onc - 
petitive in nature. Thewex.dert Cime Saving imenanecmn, fan 
plication of this technique was not evaluated. However, a 
conservative estimate is that execution time was reduced 

by a factor of five and, therefore, makes the nominal ad- 
ditional storage requirement well worth while. The specific 


operation just described is contained in SUBROUTINE REFER. 
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IV. SUMMARY AND CONCLUSIONS 


The Direction Findungesleorttnn Cansbewaomired nema 
Variety of concave or convex differentiable functions to 
find an optimal solution. However, the domain of these 
functions must be represented as a closed convex. set. 
ies applicability in practical problems ts particularly 
Meetul when the objective is te optimize alleeatren) of 
resources. 

In large scale problems it is easily adaptable to 
machine utilization. Time and stcrage requirements for a 
machine solution are primarily dependent on the nature of 
the specific problem being solved. The algorithm has 
Mmeelacively little effect onsindehine storage requirements. 
mewever, the algorithm’s impact on machine time 1S the re- 
sult of the boundary problem previously discussed. A 
method of overcoming the boundary problem as well as an 
eeetended applicatitem of this algorithm into an area of 
problems which require the allocated variable to be repre- 
Semted in typical matrix form rather than in vector torm 
1s currently being investigated by the aforementioned 


ie. Danskin. 
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APPENDIX A: FLOWCHART-GENERAL LOGICAL PROCEDURE FOLLOWED 


FIND DIRECTION 
Ol lb | 


. ENTER STARTING 
SOLUTION 


CALCULATE 2's 


PNGREAS E 

















ATTEMPT MOVE 
IN DIRECTION 
OF WASTES! 
NCIS 


ADJUST 
DISTANCE 





HAS 
OBSECTIVE 
FUNCTION 

MPROVED? 


YES 


MOVE TO NEW 
ALLOCATION 
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APPENDIX Be oUEROU TT MVE Sao ihks 


SUBROUTINE AMMCAL: AMMCAL evaluates the objective function. 
SUBROUTINE AREA: AREA determines which j's are Pee Gated 
With a given i in the ae term. 

SUBROUTINE CHANGE: CHANGE converts the optimal solution 
obtained from the 400 grid case into the starting solution 
for the 3600 grid case. It also changes values of various 
parameters used within the program for use in the 3600 grid 
solution. 

PUOBROULINE CNTROL: CNITROL controls the sequencing of events 
in solution of the problem. It tests each attempted move 
for success or failure and adjusts distance (D) accordingly. 
SUBROUTINE DF$ALL: DF$ALL is the Direction Finding Algo- 
mechm which finds the direction of fastest incréase for the 
emeare allocation vector. 

SUBROUTINE DF$LST: DF$LST is the Direction Finding Algo- 
rithm which finds the direction of fastest increase just 
for those variables on the list. 

SEMEROUTINE LIST: LIST compiles a list of those variables 
which are positive and/or those variables which have a 
positive y. 

SUBROUTINE MOVE: MOVE adjusts the allocation vector to 

mer lect di detemptred move in ee anoeeaare Oe Paso Ss tain 
crease. It moves to the boundary any variables which are 


ae 
within 10 Cre ened bo OUine anya 


Be 





SUBROUTINE OMGALC: OMGALG calculates the veetorn or fice 
pattial derivatives omeencwalloecatrons Voetor: 

SUBROUTINE REFER: REFER fills a matrix with values rep- 

resenting geometrical distance between two grid : deme 

SUBROUTINE SET$P: SET$P sets precision desired into the 


program and fills the P vector. 
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BELOW IS RIGHT SINE OF OPTIMAL SOLUTION OF 400 GRID CASE 
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IF(OMM.LT.~AQM) AOM=O0MM 
520 Cee ete 


U=AOM 
NOW TO CHECK FOR ELEMENTS ON UPPER BOUNDARY AND NOT IN DS 


540 JMARK=0 
DM 550 JJ=1,LONGS 
TECINDCS)) GO-1T0 S59 
LEGS tyes 0) GO TD 550 
IF(OMS(JSJ).G6E.A0M) GO TO 550 
 AOM=0MS (SJ) 
JIMARK= 39 
550 CONTINUE 
LP een) CETT0 570 
NOW TO RECALCULATE MU 
MU=(N*MU+F+A0NM)/ (N41) 
AOM=MU 
N=N+1] 
IND{ JMARK)J=.TRUE. 
GO TQ 549 
NOW TO CHECK FOR ELEMENTS ON LOWER BOUNDARY AND NOT IN ODS 
560 JMARK=0 
DO 565 JJ=1,LONG$ 
Pet iNet ss)) Gao to 565 | 
Wet tae). Gl.060) GO 18-565 
TFC(QMS( SS) LE-AGM) GQ TO 565 
AOM=0M5 (JJ) 
JIMARK=JJ 
565 CONTINUE 
Pagan Ls) GO TM 880 
RECALCULATE MU AGAIN 
MU=(N*XMU4+A0M)/ (N41) 
AOM=KMU 
N=N4+1 
IND ( JMAQKI=.TRUE. 
GO TM 560 
S7TO TFC(KK.LFO.1L) GI TO 590 
KK=1 IMPLIES NO MORE SEARCHENG IS REQUTRED 
GO TO 560 
580 KK=l 
GQ TA §40 
NOW TO FIND SUM OF SQUARES OF ELEMENTS IN DS 
S90 S5-0.0 
DO 595 JJ=1,LONGS 
Bei smote INOC IS) ) GO TO 595 
SO=CNS$ (SSI -—MU 
SS=SS+#+S0*SN 
595 “GONTINUE 
MEA SS.LE.VS) RETURN 1 
LA=DSORT(SS) 
LAMBDA EQUALS SQUARE RQOT FF SUM OF SQUARES 
DD 46C0O JJ=1L,LONGS 
Beier) eINOCJIJ)}) GA ¥O £00 
$(1) NOT IN DS IMPLIES GAMMA(T)=0 
SO=OMS(SS)-MU 
GAS( JJ) =SD/LA 
GAMMA VFCTOR IS FILLED HERE 
600 CONTINUE 
RETURN 
END 
SUBRDUTINE AMMCAL 
THIS PART CALCULATES THF HFS AND AMM 
TMPLICIT REAL *R(A-H,C-Z,$) 
PBGICAL SHEUSTH,FAILI 
DIMENSTON H(32600),9(24500) 
COMPON/WORKBR/4(32609),0(32600) 
COMMONS XFERIL/L ee MyLLE yl @eML,yMM 
300 00 305 T=1,LCNG 
305 H(1T)=0.0 
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